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A Lagrangian flow network is constructed for the atmospheric blocking of eastern Europe and western Russia 
in summer 2010. We compute the most probable paths followed by fluid particles which reveal the Omega- 
block skeleton of the event. A hierarchy of sets of highly probable paths is introduced to describe transport 
pathways when the most probable path alone is not representative enough. These sets of paths have the shape 
of narrow coherent tubes flowing close to the most probable one. Thus, even when the most probable path is 
not very significant in terms of its probability, it still identifies the geometry of the transport pathways. 
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Eastern Europe and Western Russia experienced 
a strong heat wave with devastating consequences 
in the summer of 2010. This was due to an atmo¬ 
spheric blocking episode that lasted during sev¬ 
eral weeks. Despite these type of events have 
been well-investigated over the years, a complete 
understanding and prediction is still missing. In 
this work we present a characterization of this 
flow pattern based on the study of fluid transport 
as a Lagrangian flow network, so that the method¬ 
ology of complex networks can be applied. In 
particular, the most probable paths linking nodes 
of this atmospheric network reveal the dominant 
pathways traced by atmospheric fluid particles. 


I. INTRODUCTION 

Lagrangian analysis of transport in fluids, in particular 
in geophysical and time-dependent contexts, has experi¬ 
enced intense developments in the last decades. These 
can be roughly classified in three classes: Some of the 
approaches search for geometric objects lines, surfaces, 
usually related to invariant manifolds - which bound 
fluid regions with different properties 1-3 . In the second 
type of approaches one computes different types of Lya¬ 
punov exponents and other stretching-like fields in the 
fluid domain 4-7 . Finally, set-oriented methods 8-12 ad¬ 
dress directly the motions of finite-size regions. 

Most of these techniques focus in identifying proper 
Lagrangian Coherent Structures 13 ^ 15 , understood as bar¬ 
riers to transport or coherent regions with small fluid 
exchange with the surrounding medium. Much less is 
known about the actual routes of transport , the domi¬ 
nant pathways along which fluid particles travel and fluid 
properties are interchanged. 

In principle, the pathways are simply given by trajec¬ 
tories starting from the desired initial conditions. This 
is true when the advection dynamics is represented by 
a deterministic dynamical system and the initial condi¬ 
tion is precisely fixed. In many applications however, 


particularly in geosciences, stochastic components are 
added to the motions to better represent unresolved spa¬ 
tial scales 16-18 . Also, imprecisely stated initial conditions 
will develop into a divergent set of possible trajectories, 
because of the inherently chaotic character of advection 
by nearly any nontrivial fluid flow, particularly when it 
is time-dependent. In fact in real experiments such as in 
the deployment of buoys or balloons the trajectories of 
closely released objects diverge soon 19-21 . The so-called 
spaghetti plots 18 provide a visual representation of this 
dispersion. But they become, when many trajectories are 
represented, cluttered and unclear. Some type of cluster¬ 
ing or the selection of relevant trajectories is needed to 
highlight which are the dominant routes among a large 
set of possible trajectories. 

We have recently developed 22 a formalism that com¬ 
putes, in unsteady flows, the optimal fluid paths starting 
at given initial conditions and also optimal paths con¬ 
necting pairs of points. By optimal we refer to the paths 
which are more likely to be followed, in a well-defined 
sense made explicit below, by the fluid particles initial¬ 
ized in a finite neighborhood of the initial locations. By 
this reason they are called most probable paths. The 
methodology builds on the set-oriented techniques 8-12 
which discretize space to provide a coarse-grained de¬ 
scription of transport, and draws analogies with network 
theory 22-26 , for which tools to compute optimal paths in 
graphs are well developed. A related formalism address¬ 
ing optimal paths in time-independent flows in continu¬ 
ous time has been developed by Metzner et al. 27 . The 
optimal paths provide the main pathways or skeleton of 
the transport process in a given geographical area. Be¬ 
cause of the implicit stochastic ingredient in the coarse- 
graining procedure of set-oriented methods, this method¬ 
ology, at variance with other ones more tied to the theory 
of smooth dynamical systems, can be applied equally well 
to cases of deterministic transport and to strongly diffu¬ 
sive situations. 

In this paper we compute optimal transport paths for 
the atmospheric circulation during a blocking event oc¬ 
curring in Summer 2010 (in particular we focus our study 
for the period 20th July - 30th July) over Eastern Eu- 
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rope and Russia. This atmospheric flow has very different 
temporal and spatial scales, and is much more diffusive, 
than the oceanic flow analyzed in Ser-Giacomi at al. 22 . 
We give a more detailed description of the methodology 
sketched in that reference, and generalize it to extend 
the concept of most probable path to a hierarchy of sets 
of paths characterized by an increasing probability. The 
spatial coherence of these sets is also discussed. 

The paper is organized as follows: In Sect. II we 
summarize the definition and construction of the opti¬ 
mal pathways as most probable paths in a flow network. 
In Sect. Ill we extend this concept to sets of highly 
probable paths and give rules to establish their signifi¬ 
cance and spatial coherence. Sect. IV describes the at¬ 
mospheric blocking event, the data and models we use to 
compute the Lagrangian trajectories, and construct the 
flow network from them. Sect. V contains our results: 
optimal pathways for different dates and locations, and 
also a discussion of the statistical representativeness of 
the optimal paths on the sets of highly probable paths. 
The final Section summarizes our Conclusions. An Ap¬ 
pendix applies our formalism to a simple model flow, an 
analytic double-gyre system, so that the properties of the 
optimal and highly probable paths computed for the at¬ 
mospheric dynamics could be more easily understood in 
this simplified framework. 


II. OPTIMAL PATHS FROM LAGRANGIAN FLOW 
NETWORKS 


stance, using time-ordered graphs 22,28 . 

A fundamental assumption we make is that of a Marko¬ 
vian dynamics, i.e. at each time interval the ideal fluid 
particles are initialized with uniform density in each box, 
thus without keeping track of the trajectories at the pre¬ 
vious time step. The effect of such assumption is to in¬ 
troduce diffusive effects in the dynamics even when the 
original equations of motion are fully deterministic 29 . In 
the limit of very small boxes and very short time steps, 
this computational diffusion is suppressed and we ap¬ 
proach the perfect Lagrangian motion under the given 
velocity field (which itself can contain diffusive or fluctu¬ 
ating terms). 

In our network approach spatio-temporal particle tra¬ 
jectories are mapped into discretized paths between the 
network nodes. We define an M-step path ji between 
nodes / and J as the ordered sequence of (M + 1) nodes, 
p = {I, k \,..., kM- 1 , J}, crossed to reach node J at time 
t-M starting from node I at time to- Under the Marko¬ 
vian hypothesis we can associate a probability to each of 
these paths as 
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Our approach to find optimal paths in time-dependent 
fluid flows first represents the fluid transport dynamics 
as a time-dependent flow network 26 and then uses graph- 
theory techniques to extract from it these optimal paths. 
Following the set-oriented methodology 8-12,26 we proceed 
first by a discretization of the spatial domain of interest, 
dividing it into N non-overlapping boxes. In terms of 
the network-theory approach to transport 22,23,26 each of 
these boxes will represent a single network node. A large 
number of ideal fluid particles is released in each box. 
Under advection by a given velocity field, links between 
nodes are established by studying the Lagrangian trajec¬ 
tories of the particles exchanged among each pair of net¬ 
work nodes. This is conveniently done with a temporal 
discretization, i.e. we consider the dynamics restricted 
to a time interval [tojtjw] and divide it in time steps of 
length r, ti = t 0 + It , l = 0,1 For each time 

interval [ti-i,t{\ we integrate the equations of motion of 
each ideal fluid particle and keep track of each trajec¬ 
tory. The transport dynamics will then be described by 
adjacency matrices A (d , (l = in which a matrix 

element is given by the number of particles initial¬ 
ized at time f/_i in node I that end up at time ti in node 
J. Since the velocity field will vary in time the adjacency 
matrices will depend on the time interval considered. The 
weighted network we build will therefore have an explicit 
time-dependent character and can be analyzed, for in¬ 


is the probability of a fluid particle to reach node ki at 
time ti if it was initialized at time ti-i in node 
estimated as the ratio of the number of particles doing 
so to the total number of particles released at the initial 
node and time. The quantity Soltik) = Ej A kj is called 
out-strength of node k during the Z-th time step. 

Among all possible M-step paths between node I and 
J the one associated with the highest probability in 
Eq. (1) is called the most probable path (MPP) and 
is denoted by r/fj. Since this path depends explicitly 
on the number M of steps considered, it could be also 
named “fixed-time most probable path”. Its probability 
is denoted by P/| = ma x-^KPij)»}■ T° find the MPP 
and its probability we use an adaptation of the Dijk- 
stra algorithm 30 which takes into account the layered 
and directed structure of our time-ordered flow graph. 
The simplest implementation of the algorithm would in¬ 
volve finding maxima by searching over the full network, 
which can be a computationally expensive task. This 
is greatly facilitated by using the concepts of accessibil¬ 
ity and accessibility matrices 31 . Thus, for given / and 
J, our implementation of the algorithm consists of two 
main parts. In the first part one builds the tables Uof 
nodes accessible from I and J at time step l, i.e. the set 
of nodes which can be crossed at t = t; coming from I 
and proceeding towards J (see Fig. 1). Technically, this 
is done by including in U d ) the nodes ki for which the 



3 


two following conditions are satisfied: 
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In the second part of the algorithm one recognizes that 
the structure of expression (1) allows to maximize it by 
recursively maximizing over k\,k 2 , •••, kM— l- This is done 
by finding, for each accessible node fc; € (and only 
for them, without the need of scanning the remaining 
nodes in the full network), the highest probability P l Iki 
of the path connecting I and ki and the actual path as¬ 
sociated. For l = 1, i.e. for the first time step, trivially 
we have Pj ki = For l = 2,3, — 1 we apply 

recursively the formula 

= n ^ x (-^/fciT (I+1) fe i fc i+1 ) • (4) 


until the final point k>M = J is reached, and the maxi¬ 
mum probability, together with the associated path, are 
obtained (See Fig. 1). The same procedure can then be 
applied to any other pair of nodes (/', J'). 

Raising the number M of steps we observe a fast 
increase in the number of paths connecting two given 
nodes. It is thus crucial to understand how much the 
MPP is representative of the large set of possible paths 
joining two nodes. To assess in a quantitative way this 
issue we introduce the following quantity 
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which determines the fraction of probability carried by 
the MPP with respect to the sum of probabilities of all 
paths connecting nodes / and J. Note that the denom¬ 
inator can be simply computed as the matrix-product 
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FIG. 1. Schematics of the algorithm to find the MPP of M 
steps between I and J. a) First part: determination of the 
accessible nodes. Point A is reachable from I at t = t 2 but 
it is not possible to reach J from it in the rest of the time 
interval. Point C is not reachable from / at t = tju-i even 
if J can be reached from it. Point B satisfies both accessi¬ 
bility conditions, therefore, in contrast to points A and C, it 
belongs to the accessibility set and it will be considered in 
the calculation of the MPP. Systematic identification of all 
accessible nodes is done my applying the criteria in Eq. (3). 
The rest of the figure illustrates the recursive maximization 
procedure given by Eq. (4): b) In the first time step one as¬ 
signs to the links towards the nodes Ai and A 2 (considered 
to be the only ones in the accessibility set Ujj) the prob¬ 
abilities and T ( 1 ) ja 2 , respectively, c) For node B 1 

one considers the links from Ai and A 2 , evaluates the path’s 
probabilities T^'/AiT'^AiB! and T (1 ) ia 2 T (2 ^a 2 Bi , and se¬ 
lects the maximum one (in the figure the corresponding to 
the path I,AzB \, red lines). One repeats this for all nodes 
Bi,B 2 ,Bs in the accessibility set Uto obtain the MPPs 
between I and these nodes, and then the procedure can be 
iterated again for the accessible nodes at time tz. 


III. SETS OF HIGHLY PROBABLE PATHS 

For large values of M. the MPP progressively loses 
dominance and, on average, does not carry a significantly 
high fraction of probability. However the dynamics, char¬ 
acterized by a high number of paths connecting initial 
and final points, can be still described by a few of them, 
which together have a non-negligible probability. To see 
this we can relax the definition of MPP and define a fam¬ 
ily of subsets of highly probable paths (HPP) holding 
most of the probability. In our formulation each subset 
ICjj(r, e) is characterized by a rank 0 < r < M — 1 and 
a threshold parameter 0 < e < 1. Ideally the sets would 
contain all the paths whose probability is larger than 
ePj j. But since exhaustive searching of all such paths be¬ 
comes computationally prohibitive except for very small 
M, the second parameter r is introduced to determine 
the number of constraints imposed in the search for these 


relevant paths. Given the initial (J) and final (J) points 
we fix r nodes at intermediate times and look for paths 
between I and J made of segments which are MPP’s con¬ 
necting these intermediate nodes, by using the algorithm 
above. Different locations and times for these r inter¬ 
mediate nodes are scanned and paths with probability 
larger than eP^j are retained and incorporated into the 
set K,fj(r,e). For e —> 1, independently on the rank (or 
for r = 0) only the MPP is retained. = M — 1, e) 

contains all the paths with probability larger than ePfj. 
However, evaluation of these sets of HPPs can be com¬ 
putationally costly for high values of r, since the algo¬ 
rithm scales exponentially with r. Nevertheless interest¬ 
ing results can be obtained considering already low-order 
HPPs, i.e. r = 1 and r = 2. 

Once one of the subsets is computed we can estab¬ 
lish its significance by defining an extension of expression 
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Eq. (5): 
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where the sum in the numerator is over all the paths in 
the subset JCfj(r,e) and the one in the denominator is 
over all paths connecting I to J. 

Another important aspect of the sets of HPPs is to 
establish how close, spatially, are they with respect to 
the corresponding MPP. This is obtained with an aver¬ 
age distance function. Given two generic paths between 
initial and final points / and J, n\ = and 

H 2 = {I,h •••) J} we define their average distance as 


1 M—l 

d(/i 1,M2) = M _ 1 
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where d(ki,li) is a metric determining the distance be¬ 
tween two given nodes of the network. For a geophysi¬ 
cal transport network the geographical distance (on the 
sphere) between the centers of the nodes is the most nat¬ 
ural choice. For a given pair of nodes (J, J) the average 
distance between the subset ICfj(r , e) and the MPP con¬ 
necting them in M time steps is defined as 

V iJ = JjM ( 8 ) 

IJ /x 


where Nfj is the number of paths fi in the subset 
K.jj(r,e), and the sum is extended over all paths in the 
subset (remember that r/fj denotes the MPP). This quan¬ 
tity provides an estimation of how much paths in the 
subset deviate spatially from the correspondent MPP. A 
large deviation means that the probability to reach J 
from I is spatially spread in a large region and indicates 
furthermore the importance of considering the HPP sub¬ 
set instead of only the MPP. Small values of Vj^ imply 
HPP sets with the shape of coherent narrow tubes around 
the MPP, so that the MPP already characterizes the spa¬ 
tial pathways, even if its probability is not large. 

In the next Sections we apply the above formalism to 
the atmospheric flow occurring over Eastern Europe in 
Summer 2010. Computations of optimal paths and their 
sets in an analytic double-gyre system, a much simpler 
flow in which path properties could be more easily ap¬ 
preciated, are contained in the Appendix. 


ST & Z500 24/07/2010 12h 



FIG. 2. Geopotential height at 500hPa (contours, in m) and 
temperature (color code, in degrees C) over the region of in¬ 
terest, on July 24th, 12:00 UTC. 


A. Event description 


Eastern Europe and Western Russia experienced a 
strong, unpredicted, heat wave during the summer of 
2010. Extreme temperatures resulted in over 50000 
deaths and inflicting large economic losses to Russia. The 
heat wave was due to a strong atmospheric blocking that 
persisted over the Euro-Russian region from late June to 
early August 32 . During July the daily temperatures were 
near or above record levels and the event covered Western 
Russia, Belarus, Ukraine, and the Baltic nations. Physi¬ 
cally, the origins of this heat wave were in a atmospheric 
block episode that produced anomalously stable anticy- 
clonic conditions, redirecting the trajectories of migrat¬ 
ing cyclones. Atmospheric blocks can remain in place 
for several days (sometimes even weeks) and are of large 
scale (typically larger than 2000 km). In particular, the 
Russian block of summer 2010 was morphologically of 
the type known as Omega block that consists in a combi¬ 
nation of low-high-low pressure fields with geopotential 
lines resembling the Greek letter SI (see Fig.2). Omega 
blocks bring warmer and drier conditions to the areas 
that they impact and colder, wetter conditions in the up¬ 
stream and downstream 33 . We study the concrete period 
extended from the July the 20th to July 30th. 


B. Data 


IV. A NETWORK OF ATMOSPHERIC FLOW OVER 
EASTERN EUROPE IN SUMMER 2010 

In this section we describe the physical characteristics 
of the atmospheric event, the data used and the model 
we employ to obtain the air particle trajectories. 


Atmospheric data were provided by the National Cen¬ 
ters for Environmental Prediction (NCEP) Climate Fore¬ 
cast System Reanalysis (CFSR) through the Global Fore¬ 
cast System (GFS) 34 . This reanalysis was initially com¬ 
pleted over the 31 year period from 1979 to 2009 and 
extended to March 2011. Data can be obtained with a 
temporal resolution of 1 hour and a spatial horizontal 
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resolution of 0.5° x 0.5°. The spatial coverage contains 
a range of longitudes of 0°E to 359.5°E and latitudes of 
90°5 to 90° N. 

The variables needed as input to the Lagrangian dis¬ 
persion model described in the next section include dew 
point temperature, geopotential height, land cover, plan¬ 
etary boundary layer height, pressure and pressure re¬ 
duced to mean sea level, relative humidity, temperature, 
zonal and meridional component of the wind, vertical ve¬ 
locity and water equivalent to accumulated snow depth. 
All these fields are provided by CFSR data on 26 pressure 
levels. 


C. Lagragian Particle Dispersion Model FLEXPART 



As mentioned, the idea is to obtain the effective veloc¬ 
ity field felt by any fluid particle. Then the Lagrangian 
dispersion model (see next subsection) will integrate it to 
provide as output the three-dimensional positions of the 
particle at every time step. 

The numerical model used to integrate particle veloci¬ 
ties and obtain trajectories is the Lagrangian particle dis¬ 
persion model FLEXPART version 8.2 35 ’ 36 . FLEXPART 
simulates the long-range and mesoscale transport, diffu¬ 
sion, dry and wet deposition, and radioactive decay of 
tracers released from point, line, area or volume sources. 
It most commonly uses meteorological input fields from 
the numerical weather prediction model of the European 
Centre for Medium-Range Weather Forecasts (ECMWF) 
as well as the Global Forescast System (GFS) from NCEP 
(the one used in our study). Trajectories are produced 
by integrating the equation (the input velocity data are 
interpolated on the present particle position): 

^=v(X(i)), (9) 

with t being time, X the vector position of the air par¬ 
ticle, and v = v + v 4 + v m is the wind vector. FLEX¬ 
PART takes the grid scale wind v from the CFSR, but 
complements it with stochastic components v 4 and v m 
to better simulate the unresolved turbulent processes oc¬ 
curring at small scales. The turbulent wind fluctuations 
v* are parametrized by assuming a Markov process via a 
Langevin equation, and the mesoscale wind fluctuations 
v m are implemented also via an independent Langevin 
equation by assuming that the variance of the wind at 
the grid scale provides information on the subgrid vari¬ 
ance. Variables entering the parametrizations are ob¬ 
tained from the meteorological CFSR fields. For addi¬ 
tional details we refer to Stohl et al. 35,36 . 


D. Network construction 

We focus our analysis on the domain in between 0°E - 
80°E and 40°N - 70°N. In order to define the nodes of the 
network we discretize this region in 626 equal-area boxes 


FIG. 3. The geographical domain considered and the dis¬ 
cretization grid defining the nodes of our flow network. 


using a sinusoidal projection. The latitudinal extension 
of each node-box is 1.5°, the longitudinal one varies de¬ 
pending on the latitude (see Fig. 3). The area of each box 
is 27722 km 2 , so that the typical horizontal size is of the 
order of 166.5 km. This is a moderate coarse-graining 
of the resolution (0.5° x 0.5°) of the NCEP data used 
for particle integration. We take r = 12 hours as time 
discretization, which is enough to follow the dynamics of 
the blocking event. It has been shown in an oceanic flow 
network 22 that the value of r has a minor influence on op¬ 
timal paths, being more important the total time-interval 
considered Mr. We uniformly fill each node with 800 
ideal fluid particles releasing them at 5000 m of height, 
a representative level in the middle troposphere. FLEX¬ 
PART trajectories are fully threedimensional, but by ini¬ 
tializing at each time-step particles in a single layer we 
are effectively neglecting the vertical dispersion (which 
is of the order of 800 m in the r = 12 h time step) and 
focussing on the pathways of large scale horizontal trans¬ 
port. Fully three-dimensional flow networks will be the 
subject of future work. 


V. RESULTS 
A. Optimal paths 

Equipped with the tools developed above we can now 
compute pathways of transport during the atmospheric 
event described in Sect. IV. Figure 4a shows all the opti¬ 
mal paths leaving a node in the Scandinavian Peninsula 
at July 25 and arriving to all nodes which are reached 
in M = 9 steps (i.e. 4.5 days). The graphical repre¬ 
sentation joins with maximal arcs the center of the grid 
boxes identified as pertaining to the MPP. The actual 
particle trajectories between two consecutive boxes are 
not necessarily such arcs. The paths are colored accord¬ 
ing to their probability value Pj j. The MPPs with high- 
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FIG. 4. Paths of M = 9 steps of r = 12 hours in our flow 
network with starting date July 25tli 2010 (panel a)) and 
July 20th 2010 (panel b)), represented as straight segments 
(in fact, maximal arcs on the Earth sphere) joining the path 
nodes. MPPs originating from a single node (black circle) and 
ending in all accessible nodes. Color gives the Pjj value of the 
paths in a normalized log-scale between the minimum value 
(deep blue) and the maximum (dark red). Panel a): proba¬ 
bilities ranging from 1CF 3 to 10~ 14 . Panel b): probabilities 
ranging from 10 -3 to 10 -15 . 


est probability (reddish colors) follow a dominant anticy- 
clonic (i.e. clockwise) route bordering the high pressure 
region (see Fig. 2, but note that this is at a particular 
time, whereas the trajectory plots span a range of dates 
of more than four days) without penetrating it. There is 
also a branch of MPPs with much smaller probabilities 
(yellow and bluish colors) that are entrained southward 
by a cyclonic circulation. 

Despite the persistent character of the Eulerian block 
configuration, sets of Lagrangian trajectories become 
highly variable in time. See for example the set of MPPs 
starting from the same initial location but five days ear¬ 
lier (Fig. 4b). The southward cyclonic branch is now ab¬ 
sent, all MPPs following initially the anticyclonic gyre. 
Remarkably, the set of trajectories bifurcates into two 
branches when approaching what seems to be a strong 
hyperbolic structure close to 40°N 60°E. A hint of the 
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FIG. 5. Optimal paths of 9 steps of r = 12 hours with starting 
date July 20th 2010, entrained in the high- and in the two 
low-pressure areas of the blocking. Same coloring scheme as 
in Fig. 4. Panel a): probabilities ranging from 10~ 3 to 10 -16 . 
Panel b): probabilities ranging from 10~ 2 to 1CF 1G . Panel c): 
probabilities ranging from 10 -3 to 10 -13 . 


presence of second hyperbolic structure is visible at the 
end of the westward branch, close to 50°N 30°E. Figure 
5 displays additional MPPs starting also at July 20th, 
but initialized inside the main anticyclonic region of the 
blocking, and in two low-pressure regions flanking it. Fig. 
5a clearly shows the main anticyclonic circulation, high¬ 
lighting also the escape routes from the high-pressure 
zone, associated with the hyperbolic regions described 
above. The other two panels show the cyclonic circula¬ 
tions at each side of the high, in a characteristic Omega¬ 
blocking configuration. It is remarkable the compactness 
of the trajectories inside the eastern low-pressure area, 
which form a very localized and coherent set with prac¬ 
tically no escape in the 4.5 days time-interval displayed. 

We stress that the plots in Figs. 4 and 5 are different 
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FIG. 6. Ranking plot in which the Pfj values of all MPPs ob¬ 
tained for M = 6,8, and 10 starting on July 25th in the whole 
area are plotted in decreasing order. The range of probability 
values of the MPPs can be read from the vertical axis (from 
a few percent to 10 -16 for M = 6 or to less than 10 -20 for 
M = 10). The total number of optimal paths can also be 
read-off from the horizontal axis. 


from spaghetti plots for which many available trajectories 
are plotted from different or related initial conditions. 
For our set of particles this will give 800 trajectories em¬ 
anating from each box. Here we are plotting just one 
path, the MPP, for each initial and final box pair, which 
strongly limits the number of paths from each box but, 
as we will see more thoroughly, it is still representative 
of the trajectories of many released particles. 


B. Relevance of the MPPs 

The range of colors in Figs. 4 and 5 indicates that, 
given an initial box, not all MPPs leading to different 
locations are equally probable. This is quantified by the 
probability Pfj which gives a weight to each MPP. In¬ 
dead Pjj takes a very large range of values. Figure 6 
shows a ranking plot in which the values of all MPPs of 
a given M and started at a particular date are plotted 
in decreasing order. We see a huge spread on the values 
of Pjj. Very low probability values arise because of the 
exponential explosion of the number of paths between 
two nodes with increasing M. Given these low values of 
Pfj except for the smallest values of M, one should ask 
how representative are the MPPs for the full set of paths. 
Figure 7a shows distributions of the parameter A fj(r,e) 
giving the relative importance of the different types of 
paths. We see that A-values are small when considering 
only the MPPs (r = 0), but the distributions shift to¬ 
wards higher values for paths sets of increasing r. Figure 
7b gives mean values of the A distributions. They de¬ 
crease with M, reflecting the lack of representativeness 
of the smallest sets of paths for large M. However, al¬ 
ready for r = 1 the set of HPPs has a mean value higher 
than 0.5 for a relevant range of time steps. 

Thus, for the values of M and e discussed here, the set 
of HPPs with r = 1 seems to be rich enough to represent 
the transport pathways. But how different is the geom¬ 
etry of the different paths in this HPP set? And how 
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FIG. 7. a) Normalized probability density /(A) of the merit 
figure A/j(r, e) of paths started on July 20th 2010 for M = 
9 and t = 0.1, with r = 0 (only the MPPs, black curve), 
r = 1 (blue) and r = 2 (green), b) Mean value of the A jj(r, e) 
distributions (paths’ starting date July 25th) as a function 
of the number of time steps M for r = 0 (only MPPs, blue 
squares), r = 1 (red circles) and r — 2 (single black star). 


different is it from the MPPs? We plot in Fig. 8 exam¬ 
ples of all HPPs with r = 1 and e = 0 for particular (/, J) 
values and dates. In all the cases the sets remain coher¬ 
ent and narrow tubes of trajectories defining roughly the 
same pathway as the MPP. 

A quantification of the width of the tubes can be done 
with the distance measure Dfj in Eq. (8). An average 
of it over pairs of locations is shown in Fig. 9. Although 
the tube width increases with M, it remains always be¬ 
low the typical linear box size of approximately 166.5 
km (see Sect. IV D) indicating that the tubes remain 
narrow. Thus we conclude that, despite the decreasing 
probability of the MPPs for increasing M, they remain 
good indicators of the dominant pathways in the trans¬ 
port network. 

As a final description of properties of the dominant 
transport paths, we present in Figure 10 (compare with 
Fig. 8d) an example on how the MPP and the HPPs be¬ 
tween a fixed pair of nodes change when considering dif¬ 
ferent values of M , defining the temporal interval. Typ¬ 
ically, the probability of the MPP shows a maximum at 
some intermediate value of M in between shorter values 
of M for which very few particles connect the two nodes, 
and larger values of M for which the increasing number 









FIG. 8. All paths in JCjj(r = l,e = 0.1) for different /, J, initial point I marked by a circle and final point J marked by a 
square. The color bar gives in logarithmic scale values ranging from the maximum one (dark red), corresponding to the 
MPP, to the minimum of O.IPjj- Panel a): M = 8 steps, with starting date July 25th 2011; Pjj = 7.8 x 10 -5 . Panel b): 
M = 12 steps, with starting date July 25th 2011; Pjj = 2.7 x 10 -5 . Panel c): M = 11 steps, with starting date July 20th 
2011; Pjj = 7.4 x 10 -7 . Panel d): M = 11 steps, with starting date July 20th 2011; Pjj = 1.5 x 10~ 7 . 
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FIG. 9. Plot of the mean distance Vjj (Eq. (8)) as a function 
of M for r = 1 and e: = 0.1. The quantity is further averaged 
over all the HPPs starting on July 25th. Units are kilometers. 


of factors smaller than one in the product (1) defining 
(pjj)n makes this quantity to decrease again until van¬ 
ishing. For the example shown in Figs. 8d and 10, the 
value of M giving the maximum P^j is around M « 9, 
i.e. Mt = 4.5 days. Note that the HPP trajectories 
change length but keep a similar shape in the range of 
M considered, indicating that in this time interval the 
blocking atmospheric structures evolve slowly. 


VI. CONCLUSIONS 

We have introduced MPPs and sets of HPPs as tools to 
visualize and analyze dominant pathways in geophysical 
flows. We have computed them for an atmospheric block¬ 
ing event involving eastern Europe and Western Russia. 
The computed optimal paths give a Lagrangian view of 

















9 



1.0 

0.9 

0.8 

0.7 

0.6 

0.5 

0.4 

0.3 

0.2 

0.1 

0.0 



1.0 

0.9 

0.8 

0.7 

0.6 

0.5 

0.4 

0.3 

0.2 

0.1 

0.0 


FIG. 10. All paths in K.jj(r = l,e = 0.1) for the same I,J 
and starting date as in Fig. 8d. Same coloring scheme as in 
Fig.8. Panel a) M = 7 steps; P/j = 2.2 x 10" 5 . b) M = 9 
steps; P" = 2.3 x 1(T 4 . 


the Omega-block configuration, with a central anti- cy¬ 
clonic circulation flanked by two cyclonic ones. Moreover 
they give additional insight on it, such as the variability 
of the dominant pathways, and the identification of es¬ 
caping and trapping regions. The statistical significance 
of single MPPs decreases with the time interval consid¬ 
ered, but we find always that the MPPs remain repre¬ 
sentative of the spatial geometry of the pathways, in the 
sense that the sets of HPPs are coherent narrow tubes 
providing transport paths always close to the optimal 
path. This spatial coherence of transport between pairs 
of locations was already noticed in an ocean flow 22 and 
it is also present in the model flow discussed in the Ap¬ 
pendix. Then, it seems to be a general characteristic of 
flow networks. 
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APPENDIX A: OPTIMAL PATHS IN A SIMPLE MODEL 
SYSTEM 

In this Appendix we display optimal paths and sets 
of optimal paths for an analytic model flow, the double¬ 
gyre. See for example 37,38 for basic properties of this sys¬ 
tem and computations of its Lagrangian coherent struc¬ 
tures and Lyapunov fields. Because of the simplicity 
of this flow as compared with the atmospheric situation 
studied in the main text, characteristics of the optimal 
paths could be appreciated more easily. 

The double-gyre is a two-dimensional time-periodic 
flow defined in the rectangular region of the plane x = 
{x,y) € [0,2] x [0,1]. It is described by the streamfunc- 
tion 

*p{x, V,t)=A sin( 7 r f(x, t)) sin(Try) , (Al) 

with 

f(x,t)=a(t)x 2 + b(t)x (A2) 

a(t) = 7 sin(o;t) , (A3) 

b(t) = 1 — 27sin(wf) . (A4) 

From these expressions, the velocity field is 

dib 

x = ——— = — 7 rAsin( 7 r/(a:,t)) costny) (A5) 

dy 

dib ... . .df(x,t) 

y = — = nA cos{nf(x,t)) sm(ny )^-^— . (A 6 ) 

For 7 = 0, this flow is steady. Ideal fluid particles fol¬ 
low very simple trajectories: they rotate following closed 
streamlines, clockwise in the left half of the rectangle, and 
counterclockwise in the right one. The central streamline 
x = 1 , a heteroclinic connection between the hyperbolic 
point at ( 1 , 1 ) and the one at ( 1 , 0 ), acts as a separatrix 
between the two regions. When 7 > 0, more complex 
behavior including chaotic trajectories arises. The pe¬ 
riodic perturbation breaks the separatrix, so that now 
some interchange of fluid is possible between the left and 
the right part of the rectangle. The geometric structures 
involved in this interchange have been studied with a 
variety of techniques 37,38 but the framework of optimal 
paths developed in this paper seems quite natural for this 
purpose. 

We take the parameters A = 0.1 and w = 2n/5, 
and compute paths in our network framework for two 
qualitatively different situations, namely the steady case 
7 = 0 , and the periodically perturbed case (of period 
27 t/lo = 5) with 7 = 0.3. We discretize the fluid domain 
into 100 x 50 = 5000 square boxes, defining the nodes in 
our flow network, and compute the adjacency matrices 
A®, l = 1...M, by releasing 400 particles from each of 
the boxes. In all the cases shown below we compute paths 
of M = 6 steps of duration r = 1, starting at t 0 = 0. 

Figure Al considers the steady flow (7 = 0) and shows 
all optimal paths emanating from two particular initial 
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FIG. Al. Double-gyre 6 -step paths in the steady case 7 = 0, 
from two different starting nodes (the black circles) to all the 
accessible destinations. The network nodes pertaining to each 
path are joined by straight line segments colored according to 
the path probability. The color code is logarithmic in the 
full probability range, which is [0.054, 7.73 x 10 -10 ] for the 
left node (there is a total of 94 paths emanating from it) and 
[0.0234,10 -7 ] for the right one (87 paths). 


FIG. A2. 6 -step paths for the periodically perturbed double 
gyre at 7 = 0.3, from three different starting nodes (black 
circles) to all the accessible destinations. Color coding as in 
Fig. Al, with probability ranges which are [0.0327,3.355 x 
10~ 7 ] (bottom-left node, 66 paths), [0.0245,1.879xl0~ 9 ] (top- 
left node, 108 paths), and [0.0128,1.335 x 10 -8 ] (right node, 
106 paths). 


nodes and reaching all nodes accessible from them af¬ 
ter the 6 steps. We see the general clockwise and anti¬ 
clockwise circulations at each side of the separatrix. The 
two halves of the domain remain isolated. Note that the 
paths are different from the closed streamlines. This is 
so because the discretization of the fluid domain into fi¬ 
nite boxes, together with the Markov assumption, intro¬ 
duces an stochastic component equivalent to an effective 
diffusivity 29 and leads to dispersion of the particles start¬ 
ing from a single node. In our atmospheric velocity flow 
there were in addition explicit stochastic terms modeling 
turbulent diffusion and mesoscale fluctuations. Note also 
that, as in the atmospheric case, a huge range of values 
of Pj'j is present. 

Figure A2 shows optimal paths for the periodically 
perturbed flow (7 = 0.3). The general clockwise and 
counterclockwise rotations still remain, but now there 
are pathways connecting the two halves of the domain. 
Note the strong divergence of close pathways when they 
approach the hyperbolic region at the bottom of the do¬ 
main, and how is this geometric structure what allows 
transport of fluid between the two regions that were iso¬ 
lated in the steady case. 

In Fig. A3 we display sets of HPPs between three pairs 
of nodes at 7 = 0.3. More specifically we compute the 
paths obtained with r = 1 and a probability larger than 
5% of the P/j for these pairs of nodes (i.e. the paths 
in the set A 7/ (r = l,e = 0.05)). The HPPs arrange in 
very narrow tubes around the MPP, which is the same 
behavior observed in the atmospheric paths and also in 
ocean calculations 22 . The central path in Fig. A3 clearly 
identifies the pathway followed by particles to connect 



FIG. A3. Three sets of HPPs in the double gyre for 7 = 
0.3. We show all HPPs in Affj(r = 1, e = 0.05), with 
starting nodes I at the black circles and destinations J at 
the black squares. Color coding as in Figs. Al and A2, 
with probability ranges [0.0282,0.00208] (bottom-left node, 
12 paths), [0.0321,0.00179] (top-left node, 7 paths), and 
[0.0128,0.000946] (right node, 10 paths). 


the left and right regions, using the “opening” around 
the hyperbolic region at the top of the domain. 
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